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ABSTRACT: This paper presents some experimental results that indicate 
the plausibility of using non-convex variational principles to reconstruct 
piecewise smooth surfaces from sparse and noisy data. This method uses 
prior generic knowledge about the geometry of the discontinuities to 
prevent the blurring of the boundaries between continuous subregions. 

We include examples of the application of this approach to the reconstruc¬ 
tion of synthetic surfaces, and to the interpolation of disparity data front 
the stereo processing of real images. 
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1. Introduction. 


In the processing of two-dimensional signals, often arises the problem of 
reconstructing a piecewise smooth surface from noisy observations taken at sparse 
locations. In this reconstruction, it is important not only to interpolate smooth 
patches ovei uniform regions, but to locate and preserve the discontinuities that 
bound these regions, since very often they are the most important parts of the 
surface. They may represent object boundaries in vision problems (such as image 
segmentation; depth from stereo; shape from shading; structure from motion, etc.); 
geological faults in geophysical information processing, etc. 

The most successful approaches to this problem [T2] consist of, first, interpolating 
an everywhere smooth surface over the whole domain; then, applying some kind 
of discontinuity detector (followed by a thresholding operation) to try to find 

the significant boundaries, and finally, to re-interpolate smooth patches over the 
continuous subregions. 

The results that have been obtained with this technique, however, are not 
completely satisfactory. The main problem is that the task of the discontinuity detector 
is hindered by the previous smooth interpolation operation. This becomes critical 
when the observations are sparsely located, since in this case, the discontinuities 
may be smeared in the interpolation phase to such a degree that it may become 
impossible to recover them in the detection phase. 

One way around this difficulty is to perform the boundary detection and 
interpolation tasks at the same time . In the method we will present, this is done 
by generating a variational principle that includes our prior knowledge about the 
smoothness of the surface and about the geometry of the discontinuities, as well 
as the information provided by the observations. The global minimum of this 
(non-convex) energy functional is then found by a stochastic approximation 
scheme. 

This approach is based on the work of Geman and Genian [Gl]. Before 
describing it, let us formulate the problem in a more precise way. 

Imagine a region fi of the plane which is formed by a number of subregions 
separated by boundaries which are known to be piecewise smooth curves. Suppose 
that within each of these subregions, some property / (in what follows, we will refer 
to / as "depth”) varies in a smooth fashion, presenting, at the same time, abrupt 
jumps across most of the boundaries. Suppose also that we have measurements for 
the values ot / at some discrete set of sites S\ these measurements will, in general, 
be corrupted by some form of noise. 

Our problem is then to estimate the values of / on some finite lattice of points 

T C f2, and to find the position of the boundaries, using all the available information 
in an optimal way. 
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Figure 1. Sites 1, 2, 3 and 4 arc the neighborhood of site j 

2. Gcman’s Work on Bayesian Image Restoration. 


2.o \ Images as Markov Random Fields. 


The first idea on which this approach is based, is that an image formed 
by regions of constant intensity separated by piecewise smooth boundaries can be 
modeled as a sample function of a stochastic process based on the Gibbs distribution, 
that is, as a Markov Random Field (MRF) [1] (sample functions of actual MRF’s 
can be found in his paper. See also [Cl] and [G2]). This means that the conditional 
probability of a given pixel j having a particular value /y, given the values of / in 
all the remaining sites of the lattice, is identical to the conditional probability given 
the values of / in a small set of sites which we will call the neighborhood of j. 


Given a system of neighborhoods on a lattice, we define a "clique” C as a 
set of sites such that all the sites that belong to C are neighbours of each other. 
For example, on a 4-connected lattice (Fig. 1), the sites 1, 2, 3 and 4 form the 
neighborhood of site /, and the cliques are sets consisting either of single sites, or 
of two (vertically or horizontally) adjacent sites (nearest neighbours; see Fig. 2). 

It can be shown [Pi], that for a MRF with a given system of neighborhoods, 
the probability density of the images / generated by this process is of the form: 

p,u) - 4 




where Z is a normalizing constant, /? is a parameter, and the "Energy function 
U(f) is of the form: 
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Figure 2. Cliques for the 4-connccted lattice of fug. 1. 

u(f) = Tv c (f) 

c 


where C ranges over the cliques associated with the given neighborhood system, 
and the potentials Vc{f) are functions supported on them, dims, in our example of 
a 4-connected lattice, U would be: 


r(/) = ETO+ E va(/i,/ y ) 

j mCuV w 

where i, j £ N n means that i and j are nearest neighbours, and V\ and Vj are some 
functions [2]. 

In particular, if we want a Markov random process that generates piecewise 
constant surfaces, we may use potentials: 


Vi(/} = 0 




if /, = /; 

otherwise 



The parameter (3 can be interpreted as the natural temperature of the system, 
and controls the expected size of the regions formed by the process, larger regions 
being formed at low temperatures. 


2.2 Bayesian Estimation of Markov Random Fields. 
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Suppose that we have observations g that can be modeled as: 

9j = n i) ies ( 2 ) 

where ny is a white noise process uncorrelated with /; S is some set of sites; fly is an 
operator with local support (representing blurring, for example), and is invertible 
with respect to ny. (4» may represent, for example, noise addition or multiplication 
followed by a pointwise transformation). 

The conditional probability P(g \ f) is given by: 

p (<>. I f) = R PnirHtj - h,v)) 

ies 


where P„ is die probability density function of the noise. 
The posterior distribution is found by Bayes rule: 


P(f 1 1) = 


P/(f) ■ Pjg I /) 

m 


Replacing the expressions for P f (f ) and P{g | /), taking logarithms, and 
remembering that P{g) is a constant tor a given set of observations, we get that the 
Maximum a Posteriori (MAP) estimate for / is found by minimizing: 

E U) = \u(f) - Y. ln[P n ('p- 1 ( 3 y - fly(/))] (3) 

p ies 


In particular, if n is a zero-mean Gaussian white noise stationary process with 
power spectral density a 2 , and 


Qj — Hj(f) + n j — fj + n j , 


then, 


P(9 I /) = 



1 

2^ 


E fe - fj) 2 ] 

jes 


For our example of piecewise constant surfaces, with potentials given by (1), 
the MAP estimate will be obtained by minimizing: 


E (f)= E Vaifi.fi) + 

i,jeN N 


2 <T 2 


E (fj — 9if 

j€S 



As we can see, this expression has two terms: One that measures the agreement 
of the estimate with the observations, and another that corresponds to the constraints 
imposed by our prior knowledge about the nature of the solution. The tradeoff 
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between them is controlled by the parameter 



which corresponds to the signal to noise ratio. 

It is interesting to note that the general form of this expression is similar to 
the variational principles obtained by regularization methods for solving ill-posed 
problems, playing the role of the regularization parameter. In the standard 
regularization methods, however, the functional is convex (because of the choices 
of norms and stabilizing functionals, see [P2,T1]), whereas in the case of equation 
(4), the high non-linearity of V 2 makes E{f) non-con vex. For this reason, its 
minimization becomes computationally much more expensive. 


2.3 Simulated Annealing and the Minimization of E(f). 


In a recent paper, Kirkpatrick, et. al. [Kl] proposed a stochastic approximation 
method for solving combinatorial optimization problems, which can be used for our 
minimization problem. It is based on an algorithm invented by Metropolis [Ml] 
that simulates the behaviour Gf many -particle systems in thermal equilibrium: 


Consider a system with N particles, each of which may be in any one of 
a discrete number of allowable states. Let fj denote the current state of the j th 
particle, 7, the temperature, and let E(f) be the total energy of the system. Suppose 
that we visit the particles of the system in some random sequential order. When a 
particle j is visited, we update its state as follows: 


(i) Choose a new state fj randomly from the set of allowable states (excluding 
the current state fj) , using a uniformly distributed random number. 

(ii) Compute the increment in energy A E, that results from moving the state 
of the j th particle from fj to /•. 


(iii) If AEj < 0, make the move, i.e., set fj — fj. 

If AEj > 0, generate a new random number r, uniformly distributed 
between 0 and 1. 


If r < e ~ AE */ T , set fj = fj. 

If r > e~ AE >/ T , leave fj unchanged. 

It can be shown [M1,K2] that if every particle is visited infinitely often, 
this procedure will generate global states for / distributed according to Gibbs 
distribution, i.e., as the number of iterations goes to infinity, we will have: 



As T goes to 0, this distribution will tend to an impulse (or set of impulses) 
corresponding to the state (or states) of minimum energy, that is, to the value of / 
that minimizes E(f) globally. 

One serious difficulty, however, is that attaining thermal equilibrium might take 
a very long time at low temperatures. Kirkpatrick’s idea was to start at a relatively 
high temperature (where thermal equilibrium is reached very fast), and then, to 
slowly cool the system, until ’’freezing" occurs and the state stops changing. 

Geman & Geman were able to show that if the temperature is lowered at the 

rate: 


log(n + 1) ^ ' 

where n is the number of iterations, and C is a constant, this algorithm will in 
fact converge (in probability) to the set of states of minimal energy (see also [G3]). 
Also, they proved that the non-homogeneous Markov chain that corresponds to 
the annealing procedure is (asymptotically) stationary and ergodic [3], so that we 
may use time statistics to estimate the final state. We will return to this in the next 
section. 

Unfortunately, the value of the constant C for which Geman and Geman were 
able to guarantee convergence is in general very high [4] (so that the convergence 
time becomes unpractically slow), but it has been found experimentally that a 
value of C = 3(3 log 2 (where (3 is the natural temperature of the system) normally 
produces a reasonable convergence behaviour (of the order of 1000 iterations). 

The computational efficiency of Kirkpatrick’s algorithm depends on whether 
the increment in energy A Ej associated with a change in the state of the j th variable 
is easy to compute. Fortunately, the Markov property of /, and the localization 
condition on the support of the operator Hj (equation (2)) guarantee that, for our 
problem, this will be the case, so that simulated annealing is a practical, although 
expensive way of finding the minimum of (3). 

2.4 The Line Process. 

Another powerful idea in the Geman’s work is the introduction of an 
unobservable line process l into the image model. The variables associated with this 
process are located at the sites of the dual lattice of lines that connect the sites of 
the original (pixel) lattice (see Fig. 3). These variables may be binary (indicating the 
presence or absence of a boundary between two pixels), or may take more values to 
indicate the orientation of the boundary as well. In both cases, their function is to 

decouple adjacent pixels, reducing the total energy if the intensities of these pixels 
are different. 

The introduction of this process is particularly important for the discontinuity¬ 
preserving reconstruction of surfaces from sparse observations, since it allows us to 
include into the estimation problem our prior knowledge about the geometry of 
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Figure 3. Dual lattice of line elements (sites denoted by x) 
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Figure 4. Cliques for the line process 

the boundaries between regions in an explicit way. This is done by modifying the 
energy function; the new expression is: 


£(/, 0 = £ /y, J;y) + £ Vc,(l) + -„4 £(/y 


i t jEN N 
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Figure 5. Potentials for the different configurations of a line process 
where 

0, if l,j is "on" 


fjihj) —■ 


{V 2 O';, //), otherwise 


V> is defined in eq. (1); i; :l is the line element between sites i and j, and the line 
potentials V Cl have as supports cliques of size 4, such as the one shown in Fig. 
4. Every line element (except at the boundaries of the lattice) belongs to 2 such 
cliques. The values of the potentials associated with each possible configuration of 
lines within a clique must be specified, thus, for example, if we know that straight 
horizontal and vertical boundaries are likely to be present, we may use a binary 
process, and potential values as those of big. 5 (rotational invariance is assumed). If 
we want to handle more general situations (such as piecewise smooth boundaries), 
it is necessary to allow more states for the line elements, corresponding to different 
orientations, augmenting consequently the table of values for the potentials. 
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3. Extension to the Boundary Preserving Interpolation of Piecewise Smooth Surfaces. 
3.1 Continuous Intensity Values. 


Geman & Geman successfully applied their method to the restoration of 
piecewise constant images corrupted by additive white Gaussian noise, assuming 
that the grey levels corresponding to the constant subregions were either known in 
advance, or easily obtainable from the data (as peaks of a histogram, for example). 
However, if we want to apply this construction to more practical cases, we must relax 
the piecewise constant condition, so as to include tilted planes, smooth gradients, 
and in general, piecewise smooth surfaces. 

To do this, we must allow the depth variables f 3 - to take any real value. However, 
this would require, in general, to replace Metropolis algorithm by a process capable 
of handling continuous variables (such as a diffusion. See [G2]), and in this case, 
the convergence of the annealing process is uncertain (rigorous convergence results 
have not been established yet). It is also possible, of course, to discretize the set of 
allowable values for f 3 , and to use simulated annealing in a conventional way, but 
if this discretization is fine enough to guarantee sufficient precision in the results, 
the price that one would have to pay, in terms of computational efficiency, may be 
too high. 

One way around this difficulty, is to define a "mixed” annealing strategy, 
in which the continuous variables are updated in a deterministic way, while 
the stochastic updating defined by Metropolis algorithm affects only the discrete 
variables. To do this, we replace the depth potential V f in eq. (6) by one that 
guarantees that, for any given state of the line process, the resulting conditional 
energy function E(f | /) is convex, so that the use of a deterministic gradient descent 
iteration for updating the state of the depth process is justified in some sense (this 
is not a rigorous argument; we present it only to provide some motivation for the 
definition of this strategy which, at this point, is justified only by our experimental 
results). 

This can be achieved by using a quadratic potential: 


fv l a) = 


0 , 

1 ifi - fj)\ 


if kj is "on" 
otherwise 



Since the resulting conditional energy function is of the form: 

E U I 0 = Ufi - fif + A £ (fj - 91? + K (8) 

jes 

for some positive constant K, we have that E(f | /) > 0 for any / ^ 0, and by 
an appropriate change of coordinates, it can be put in the form: 



with u £ £ ,L l and A a non-negative definite matrix. 

Now, for any t £ (0,1), and any u^v, we have, 

tE{u | l) + (1 - | l) - £(*u + (1 - t)v | l) = 

— £(l — £)(u — v) t A(u — v) > 0 

so that £(/ | l) is a convex function of /, and an iteration of a gradient descent 
algorithm will move towards a global minimum of this function. Note however, that 
there may be degenerate cases in which some region Q within which there are no 
observations is isolated from the rest of the lattice by the line process. In this case, 
any solution for which 

fj — constant, j £ Q 

will be a fixed point of the gradient descent algorithm, and although all these 
solutions have the same energy, some of them are obviously more desirable than 
others. Experimentally, we have found that a good strategy to prevent the formation 
of undesirable "islands", is to use the global minimum of E(f | / = o) (which is the 
unique fixed point of the gradient descent algorithm) as the starting configuration 
for the mixed annealing process. This configuration can be interpreted physically as 
a membrane that is coupled to die data points (observations) by means of springs 
with constants equal to and is in a position that minimizes the total energy 
stored in itself, and in the springs (see |T2]). 

3.2 The Mixed Annealing Process. 

The scheme we are proposing is as follows: 

Every global iteration (that is, for every fixed temperature), all the line and 
intensity sites are visited sequentially. When a line site is visited, its state is 
updated using Metropolis algorithm. The corresponding increment in energy A Ei 
is computed as follows: 

Let Ci and C 2 be the two cliques to which the line element belongs. Let kj be 
its current state, and kj be the candidate state (if l is a binary process, kj = l - 
otherwise, it is chosen at random from the set of allowable states different from l-j). 

Let l be the current global line configuration, and / be obtained from l by replacing 
kj by kj> We have: 


where 


AE t = VcXl) ~ V Cl (l) + VcXl) ~ VcX) + 
+sgn(/,j -- fj) 2 

sgn(i) = (- 1 , 

lo. 
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if x > 0 
if x < 0 
if x — 0 




When an intensity site j is visited, its new state f*j is obtained deterministically 
the formula: 


j|C(0,l](l Sgn + Oiqjgj 

j|G(0,l](f Sgn(/ t ’y)) -j- ctqj 



where a = and qj = 0, unless there is an observation at site j (with value 
< 7 y), in which case, qj — l. 

The temperature is lowered using (5). 


3.3 Parameter Values. 


Unfortunately, we do not have, at this point, a clean way for selecting the 
values for the parameters of the system, other than a trial and error procedure. A 
few considerations about their meaning are in order: 


The parameter a controls the degree of smoothing that one applies to the data, 
and must be related to the quality of the observations. Large values of a 10) will 
force the interpolated surface to pass through the observations, while small values 
1) will smooth away all but the largest fluctuations. 


The parameters associated with the line potentials control both the shape and 
the number of boundaries that the algorithm will detect. For example, for the 
binary line process of Fig. 5, if h is the height of the smallest jump that we want 
to consider as a boundary, and d is the largest gradient in the smooth regions, we 
must have for Vj (the potential of a straight boundary) 


d 2 <Vi< h 2 (11) 

The ratio of Vj to the potentials for the rest of the configurations, represents 
our prior knowledge about the relative likelihood of corners, T junctions, etc. 


4. Experiments 


We now present some experimental results that support the use of this approach 
for surface reconstruction and image segmentation tasks. In these examples, we 
assume that we have the following prior knowledge about the nature of the surfaces 
we are trying to reconstruct: 

(i) The region under consideration can be segmented into a small number of 
subregions. 

00 Within each subregion the surface is smooth (the gradient is less than 0.5). 

(iii) The boundaries between regions are piecewise horizontal or vertical. There 
are relatively few corners. 

(iv) The average height of the discontinuities across boundaries is greater than 

0 . 8 . 
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(v) The observations are corrupted by an additive white Gaussian noise process, 
and we have some estimate of its intensity. 

This knowledge is embodied in the model for the line process, and in the 
numerical value of the parameters. For our experiments, we have used the binary 
model of Fig. 5, with the values for the potentials for the different configurations 
indicated there. 

4.1 Experimental Results. 

In the first set of experiments, we generated sparse observation points at 200 
random locations of a 30 X 30 rectangular grid. Figures 6, 7, 8 and 9 show (with 
height coded by grey level) the observations (a); the boundaries found by the 
algorithm (b); the configuration obtained by interpolation with no boundaries (c), 
and the final reconstructed surface (d), for: 

(i) A square at height 2.0 over a background at constant height = 1.0 (Fig. 

6 ). 

(ii) A triangle, with the same characteristics (Fig. 7). 

(iii) A tilted square plane (slope = 0.1) over a constant height background 
with white Gaussian added noise (a — 0.1)(Fig. 8). 

(iv) Three rectangles at different (constant) heights over a uniform background 
(Fig. 9). 

In the second set of experiments, the same algorithm was used for a boundary 
detection / image segmentation task. In this case, we have observations (corrupted 
by white Gaussian noise) at almost every point in the lattice. The original figure is 
a 10 X 10 tilted square plane of slope 0.2 located at the center of the lattice (Fig. 
10). Note that in this case, the tilted plane cuts across the uniform background, so 
that the vertical steps at both sides have opposite signs, while the horizontal steps 
change sign at the center of the figure, where, in fact, there is no discontinuity. As 
in the previous experiments, the results speak for themselves. 

Finally, we present an example of the application of this algorithm to the 
processing of real images. We use it to interpolate disparity data, obtained along the 
zero-crossing contours of the convolution of a stereo pair of aerial photographs with 
a "Difference of Gaussians" operator, by Grimson’s implementation of the Marr- 
Poggio stereo algorithm [G4,M2] (the data array was rotated, so that assumption 
(iii) held). The results are shown in figure 11. We believe that they will improve 
when we implement the extensions to this method outlined in section 5.2. 


4.2 Convergence of the Annealing Process. 

We have found that the mixed annealing algorithm, with annealing schedule 
given by (5), eventually converges to a low energy solution; however, its convergence 
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Figure 6. (a) Observations of a square at height 2.0 over a background at height 1.0 (a white 
pixel means that the observation is absent at that point), (b) Boundaries found by the Algorithm, 
(c) Interpolation with no boundaries, (d) Reconstructed surface. 
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Figure 7. (a) Observations of a triangle at height 2.0 over a background at height 1.0. (a white 
pixel means that the observation is absent at that point), (b) Boundaries found by the Algorithm, 
(c) Interpolation with no boundaries (d) Reconstructed surface. 
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Figure 8. (a) Observations of a tilted square (slope = 0.1) over a background at height 1.0 
with added white Gaussian noise (a — 0.1) (a white pixel means that the observation is absent 
at that point), (b) Boundaries found by the Algorithm, (c) Interpolation with no boundaries, (d) 
Reconstructed surface. 
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Figure 10. (a) Observations of a tilted square (slope = 0.2) over a background at height 
1.0 with added white Gaussian noise (a — 0.1). White pixels denote missing observations, (b) 
Boundaries found by the Algorithm, (c). Interpolation with no boundaries, (d) Reconstructed 
surface. 
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Figure 11. (a) Disparity data for a stereo pair of aerial photographs.(b) Interpolation with no 
boundaries, (c) Boundaries found by the Algorithm (d) Reconstructed surface. 
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Figure 12. Snapshots of the Annealing process. 

can be greatly improved if we periodically set the state of the line process using an 
estimate of the final (lowest energy) state. 

To get this estimate, we use the crgodicity of the process ([Gl], theorem C), and 
compute time statistics about its evolution. Specifically, we estimate the marginal 
distributions of the states of every "line" cell by computing, within a time window, 
the percentage of the time for which this cell is "on". If it is more than half of the 
time, the state of this cell in the annealing network is forced to 1; otherwise, it is 
forced to 0. As the temperature decreases, die probability distribution of the global 
states becomes peaked, and the maximum probability estimate obtained from the 
marginals gets closer to the true final state. 

The results discussed in the previous section were obtained as steady states of 
this modified process, using a time window of 10 iterations (for the experiment of 
Fig. 10, an additional window of 100 iterations was used). In all cases, the steady 
state was obtained after approximately 450 iterations. Fig. 12 shows snapshots, taken 
every 50 iterations, of the state of the line process for the surface reconstruction 
problem on a 100 X 100 lattice, with 2000 observations, portraying a central square 
figure at constant height over a uniform background. 


5. Discussion. 


The results we have presented indicate the plausibility of using variational 
principles that include our prior knowledge about the geometry of the discontinuities 
of a piecewise smooth surface, to reconstruct it from sparse and noisy data, preserving 
the boundaries between continuous subregions, dhis method may also have a more 
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general significance, since most problems in early vision are ill-posed and must be 
regularized [P2]. The regularization methods used so far are based on smoothness 
assumptions, and for this reason, all have problems when dealing with functions 
that are only piecewise continuous. In this sense, the approach we have presented 
may be regarded as an extension of standard regularization methods for handling 
discontinuities (see [P2]). The algorithms presented are admittedly slow. Their 
computational efficiency, however, can be greatly improved by the use of parallel 
hardware, particularly in view of the local nature of the support for the state 
updating operations (equations (9) and (10)). It would be interesting to investigate 
the implementation of this algorithm on the "Connection Machine", currently under 
construction at the A.l. laboratory [HI] (see also [FI]). 

5.1 Relation to other Boundary Preserving Surface Reconstruction Techniques. 

The use of MRF models is closely related to the use of physical models for the 
surfaces that one wants to reconstruct. Thus, for example, Terzoupulos [ T2] proposes 
the use of a thin plate model to embody the knowledge about the smoothness 
of the surface. A threshold on the bending moment of this plate is then used to 
locate the discontinuities, and the plate is allowed to break at these points. This 
technique is computationally more efficient, since the energy functionals involved 
are convex; however the use of the method we are proposing seem to have some 
definite advantages: 

(i) From a conceptual viewpoint, it is better to perform the interpolation and 
boundary detection tasks at the same time, rather than approximating 
an everywhere smooth surface first, since this operation hides the 
discontinuities that one then tries to find in the second phase. 

(ii) In our method, the values of the parameters depend only on the average 
height of the jumps that one wants to consider as boundaries in the 
reconstructed surface , and thus, they are independent of the location of 
the observations. If these are sparsely located, the bending moment of the 
thin plate may be low, even when the discontinuity is relatively large, and 
the threshold method may fail. 

(iii) A priori knowledge about the shape, orientation and position of the 
discontinuities can be easily incorporated by choice of the potentials of 
the line process. 

(iv) The same algorithm can be used for surface interpolation, noise elimination 
(smoothing) and boundary detection. 

5.2 Extensions and Open Problems. 

This method can be easily extended to the case of piecewise smooth (not 
necessarily straight) boundaries by using a 4-valued line process to include edge 
orientation (see [Gl]). This extension will increase its practical value, and permit 
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its use in a variety of Computer Vision tasks (processing of depth, motion and 
brightness data) as well as in other fields (geophysics, medicine, etc.) 


Another interesting extension is the use of the knowledge about the position of 
the boundaries provided by some process (for example, the output of a conventional 
edge detector operating on the components of a stereo pair) to influence the 
reconstruction of a different related surface (such as the disparity surface obtained 
from a stereo matcher). This can be easily accomplished, in principle, by modifying 
the potentials of the line process, so that the presence of a related edge at a line 
location lowers the energy of the corresponding configurations. We are currently 
working on these developments. 

An important problem that has to be addressed to make this technique more 
practical, is the improvement of its computational efficiency. One needs to accelerate 
the convergence of the annealing scheme, and also, see if it is possible to exploit 
the structure of the energy functional of this particular problem, to develop more 
efficient (possibly deterministic) algorithms to find its global minimum. 


From a theoretical point of view, it would be interesting to establish in a 
rigorous way the conditions for the convergence of the "mixed" annealing strategy, 
since it may be applicable to a wider class of optimization problems. 


Another open question is the determination of the optimal values for the 
parameters of the energy functionals. It will be interesting to explore in this connection 
the use of statistical methods [B1.K3], regularization techniques p2.Tl.Wl], and 
learning algorithms [H2]. 
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Notes. 


[1] A Markov random field, a generalization of the concept of a Markov chain, is 
formally defined as follows (see [G1,K2]): 

Let 5 be a set of sites, and G = {G 5 , s £ S} be a neighborhood system for S , 
i.e., a collection of subsets of S for which: 

(0 s £ G 8 . 

(ii) s G G r if and only if r e G a . 

Let F = {F a , s e S} be any family of random variables indexed by 5, and ft 
be the set of possible configurations for F. F is a MRF with respect to G if: 
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(i) P(F = /) > 0, for all fen. 

(ii) P(F S = f s | F r = f r r s) = P(F S — f s \ F r — f r reG a ). 
for every s e S. 


[2] Note that whereas the functions defining valid conditional probabilities for a 
MRF cannot be chosen arbitrarily, and are, in general, very difficult to specify 
[B1,G1], the form of the potentials V c is not restricted in any way, and can be used 
freely to specify the required behaviour of the field f (which is what one does in 
practice). The relation between these potentials and the conditional probabilities is 
given by the following formula (which follows from Bayes rule): 


P(Fi = h | F j = ^ i ) = 


exp[-^E c:iecVc{f)\ 
Z qe Q ex Ph^ 'Ec:iec Vc{f q )] 


where Q is the set of allowable values for the state of and f q is the configuration 
which is equal to q at site z, and coincides with / everywhere else. 

[3] The annealing process can be considered as a Markov chain (the discrete time 
parameter is the number of global iterations) with states corresponding to the 
global states of the field. Since the transition probabilities (defined by Metropolis 
Algorithm) are not constant (they change with the temperature which in turn depends 
on the iteration number n), this chain is non-homogeneous. It is asymptotically 
ergodic in the sense that for any real valued function Y of the global state at time 
t, f[t)> we ha ye [Gl]: 

lim i £ Y{f(t)) = / Y( U )dP f (u,) 

where n is the set of allowable global states. This means that we can use time 
averages to estimate ensemble averages. 

[4] Geman & Genian showed that convergence can be guaranteed if C is greater 
than AT A, where N is the total number of sites in the lattice, and A is the largest 
absolute difference in energies associated with pairs of global configurations that 
differ at only one site. 
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